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The Earth's density distribution can be approximately considered piecewise continuous at the scale 
of two-flavor oscillations of neutrinos with energies about 1 MeV. This quite general presumption 
appears to be enough to analytically calculate the day-night asymmetry factor. Using the explicit 
time avergaing procedure, we show that, within the leading-order approximation, this factor is 
determined by the electron density immediately before the detector, i.e. in the Earth's crust. 
1 Within the approximation chosen, the resulting asymmetry factor does not depend either on the 

properties of the inner Earth's layers or on the substance and the dimensions of the detector. For 
beryllium neutrinos, we arrive at the asymmetry factor estimation of about —4 x 10 -4 , which is 
" at least one order of magnitude beyond the present experimental resolution, including that of the 

. Borexino experiment. 
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I. INTRODUCTION 

^: 

O-l The effect of neutrino oscillations in vacuum lies beyond the Standard Model and is thus interesting both from 
^ • the theoretical and experimental point of view. The oscillations in medium are also studied since Wolfenstein, who 
showed that the neutrinos acquire a specific flavor-dependent potential due to the coherent forward scattering on 
the matter [1]. As a result, the neutrino propagation in medium should demonstrate the conversion from one flavor 
into another (i.e. flavor oscillations), even if the vacuum mixing is negligible. This spectacular result is known as 
. the Mikheev-Smirnov- Wolfenstein effect 0] and suffices to explain the deficit of observed solar electron neutrinos [|[ . 
1 According to Mikheev and Smirnov, the leading-order estimation for the electron neutrino flux depends only on the 
t-H . points where the neutrino was born (the core of the Sun) and absorbed (the detector). However, the properties of the 
medium between these two points can also slightly affect the flavor composition of the observed neutrino flux, leading, 
in particular, to the day- night (solar neutrino flavor) asymmetry The latter effect being crucial for the entire 

flavor oscillations framework, a number of experiments were set up to catch this slight flavor composition variation 
resulting from the nighttime neutrino propagation through the Earth. In April 201 1, the Borexino collaboration did 
report the possible observation of the asymmetry at the level of I0 -3 for beryllium neutrinos, however, with a high 
uncertainty Q . Since other neutrino experiments have already yielded a considerable amount of data on the vacuum 
neutrino mixing 0, Q , the question is whether it is really worth looking forward to the experimental discovery of the 
day-night asymmetry in the nearest future. Namely, one may ask: based on the data provided by other experiments, 
how close is the theoretically predicted asymmetry to the present experimental resolution? According to the analysis 
presented below, the magnitude of the day-night asymmetry is separated from this resolution by at least one order of 



5_i _ 

magnitude. 



Within the neutrino oscillations framework, it is common to use the Schrocdinger-like equation to describe the 
spatial variations of the neutrino flavor [l|, Within the two-flavor approximation, the Schroedinger problem is 
posed in terms of the 2x2 flavor evolution matrix (operator) R(x,xq), whose elements are the neutrino flavor 
transition amplitudes after traveling from the point xq to x. The evolution equation and the initial condition read, 
respectively, 

dR(x,x ) -\Ar\nr\ ft\ 
0^ = -i\A{x)R{x,x ), (I) 

R(x ,x ) = 1, (2) 
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where the matrix Hamiltonian is a point-dependent linear combination of the Pauli matrices, 

A(x) = a(x)a 3 + b<7\, (3) 

a(x) = -cos26> + 2EV ^\ & = s in20 o . (4) 
Am 2 - 

Here, 8o is the vacuum mixing angle, E is the neutrino energy, Am 2 is the difference between squared neutrino masses, 
and V(x) = V2G F N e (x) is the Wolfenstein potential, which is proportional to the electron density in the medium N e 
and the Fermi constant Gf- The constant coefficient A = Am 2 /4E is the reciprocal neutrino oscillation length, up to 
the factor of n, 

^osc = q = T- (5) 

Am A 

Equation (JTJ) defines a one-parametric subgroup of SU{2) and hence of 50(3), the translation along x being the 
group operation. In this sense, equation (JTJ is analogous to the Bargmann-Michel-Telegdi equation Q in the spinor 
representation 10]. It is well known that the solution of the matrix linear ordinary differential equation, such as IfT]). 
can be represented as a time-ordered exponent (Dyson expansion) [111 ] . However, in the case of general N e {x) profile, 
the solution in the sense of Poincare (in terms of a number- valued series without a symbolic operation such as time 
ordering) appears to be too challenging to find. In the recent investigations, considerable progress was made in finding 
exact solutions of Eq. (|T|) in certain special cases [HI, [l3j]; nevertheless, the general approach to this kind of equation 
still remains to be approximate. 

Quite a large number of publications are devoted to the analysis of these approximate solutions. Probably the 
most effective technique for finding the approximate solutions of matrix linear differential equations, such as ([1}, is 
the so called Magnus expansion [lj] , which is a modification of the Baker-Campbell-Hausdorff formula [l5[ . This 
approach provides the solution up to any order of approximation, as well as the constraints on the remainder terms 
(la |. Unfortunately, this technique fTrl - 420j | . as well as the other general methods (see, e.g., [U|2^]), do not provide 
the way to find the solution in its explicit form, not firmly fixing the Earth model, namely, N e (x) density distribution. 
It is thus desirable to find an approximate analytical solution of equation (JXJ) , which is valid under quite general 
assumptions about the electron density profile N e (x). This idea was also developed in [23Tj25| . 

In our paper, we are not only aiming at finding the relevant approximate expressions, but also at estimating their 
accuracy and applicability domain. Namely, in section UH we find the approximate solutions for the flavor evolution 
matrix inside the Earth, and then, in section Hill we arrive at the observation probabilities for the neutrinos of different 
flavors. These probabilities are finally subjected to the averaging procedure due to the continuous observation of the 
solar neutrinos throughout the year (Sec. IIII B|) . and the results of this averaging are discussed in section ITVl The 
magnitude of the day-night asymmetry appears to be sensitive to the non-trivial structure of the Earth's crust under 
the neutrino detector, so the effect of the crust is paid special attention in section llll CI In particular, we show that our 
estimation for the neutrino day-night asymmetry is valid for beryllium neutrinos, which are observed in the Borexino 
experiment. The detailed derivation of the approximate solutions for the evolution operator is moved the Appendix 
to make the flow of the paper less complicated. 



II. THE DENSITY PROFILE AND THE EVOLUTION MATRIX 



In our investigation, we use the model electron density profile N e {x) with n — 1 narrow segments, where it changes 
steeply, separated by n wide though sloping segments. Let us call these segments cliff's and valleys, respectively. 
Let the cliffs be localized near points Xj, j = l,n — 1, namely, occupy the segments [xj,x~j] of widths Cj, where 
x^ = Xj ± £j/2. Then the valleys are [xj_-y,xj] and have the widths Lj w Xj — Xj-\. In fact, this kind of model 
is a good approximation for the Earth's density profile known in geophysics, where it corresponds to the so called 
Preliminary Reference Earth Model (PREM) pa . Inside the Earth, the oscillation lengths of neutrinos with 
E ~ lMeV are of the order of 30km, and the inequality 

€j < lose < Lj (6) 

takes place, i.e. the cliffs are narrow and the valleys are wide compared with the oscillation length. In the following, 
we will briefly call such density distribution piecewise continuous. It is worth saying here that during the night, the 
neutrino ray traverses different paths through the Earth, thus, the lengths ej and Lj vary. However, the assumption 
© holds for the most part of the night. 
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The total flavor evolution operator for such piecewise continuous density profile equals the matrix product of the 
evolution operators for all cliffs and valleys. Within each of these segments, the two small parameters arise: the first 
of them, 

r, = 2EV{ f < lCT 2 (E ~ 1 MeV) (7) 



due to the relatively small density of the Earth [2J, |28j , while the second parameter due to the piecewise continuous 
structure of the density profile, 

iosc/Lj for jth valley, , , 

€j/e osc for jth cliff. W 

The evolution matrix for each segment, as well as the total evolution matrix, can be subjected to the following 
unitary transformations (ljjj : 

R(x,x ) = Z + (x)Y(4,(x))Ro(x,x )Y' 1 (^(xo))Z-(x ), (9) 
Y(tp(x)) = cost()(x) + ia 3 sin tp(x) = exp{ia3tp(x)}, (11) 



where ui(x) = y/ a 2 (x) + b 2 is the effective oscillation wave number and ip( x ) = ^ J w(x)dx is the corresponding phase 
incursion. For the calculations which follow, it is also useful to introduce the effective mixing angle in the medium 
which is defined by the expressions 

cos20(x) = -^4, sm29(x) = (12) 

UJ(X) UJ(X) 

In terms of this angle, 

Z ± (x) = exp{±io- 2 0(x)}. (13) 

The transformation with matrices Z^(x) locally diagonalizes the Hamiltonian A(x) in the point x. It thus makes 
the complete diagonalization in the homogeneous case N e (x) = const 0. The transformation with the operator 
Y(ip(x)) isolates the effect of the medium inhomogeneity, the transformed evolution matrix i?o(a;,^o) satisfying the 
equation 

dRo ^ = -i0(x)a 2 e^^Ro(x,x o ), (14) 



where the dot over 9 denotes the gradient 



Due to the fact that the neutrino detector is homogeneous (9{x) = 0), Eq. (fT4|) has a well-defined and physically 
relevant x — > +oo limit for any fixed xq. Moreover, asymptotically convergent behavior of such systems of differential 
equations is stated by the Levinson theorem (ll| . 

In the homogeneous case, the equation above is trivial. However, in the valley [x~j~ ,xj +1 ], the slow change of the 
density N e (x) enables us to use the so-called adiabatic approximation leading to the same result 

Ro(x~ +1 , 2+) —1 + 0(?75valloy), (16) 

where the remainder term is a (generally speaking, non-diagonal) matrix. On the other hand, if the Wolfenstein 
potential undergoes a considerable change within the narrow cliff [xJ,Xj] with the phase incursion A?/> <C 2tt, then 
we get 



R 



o(x+,xj) - exp {-ia 2 A^e 2iCT ^K 7 ) j + 0( V 6 cliS ), (17) 
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where A9j = O(x^) — 9{x- ) = 0{rf) is the jump of the effective mixing angle on the jth cliff. Moreover, one can show 
that within the more accurate 0(r)6) approximation, the above expressions take the form (see the Appendix) 

Ro(xJ +li x+) = exp|-^[e 2i ^^ +1 ^(^ + 0- e 2i ^4^(x+)]| + O(^ 2 ) (valley), (18) 

R (x+,xJ) = exp{(-i<7 2 A0 3 ;+ia 1( i j )e 2ia3 ^} + 0{r,5 2 ) (cliff), (19) 

where 

't 

IH = 2A J(V ~ <cjW(y)dy = 0( V S). (20) 

x 7 

Now let us write the evolution operator for the whole neutrino path. The neutrinos observed during the day are 

created in the point Xq inside the solar core, then travel to the Earth, enter the detector in the point x± and are 
finally absorbed in the point x* inside it. In the nighttime, however, after reaching the Earth in the point Xi, the 
neutrinos pass through n — 1 Earth's layers (valleys) discussed above, and only after that they enter the detector 
in the point x n and are absorbed in x* . Crossing the Sun-to- vacuum interface, as well as traveling inside the Sun, 
does not involve steep electron density changes, thus we can treat the whole segment [xq, x±] as a single valley. As it 
was mentioned earlier, the flavor evolution operator for the whole neutrino path is a matrix product of the evolution 
operators for each segment (each valley and cliff). By the substitution of the approximate solutions ([18)) and ([19]) 
into representation @, after some transformations we find the total evolution operator in the form (compare with 
Eq. §M0§) 

R(x*,x ) = R Aet {x*^ x +y°*K e ™^-en/2\) e i^ 

x ... x e^A^g-iosAfligloi/iietosA^iguriisnn^g-iojiflsm. + 0{nr,5 2 ). (21) 

Here, the subscript 'Sun' refers to the point xo inside the solar core, where the neutrino is created, and the evolution 
operator inside the neutrino detector is denoted i?det- The factor n in the remainder term indicates that it contains 
the sum over all cliffs and valleys. Moreover, we use the following notation: 



9f = d(xf), A9j = 9(x+) - 9(xj), j = l,n-l, (22) 
9f = 9(xf), AOj = 9(x+) - 9\xJ 



■ &{xf), a9j = 6(xf) - 9(xJ), j = l,n-l, (23) 



N = H + ^ = I ( 2A(x - xj)9(x) + ^ j dx, j = 1, n - 1, (24) 



Aipj = ip(Xj ) — il>(xj_i) = A / u)(x)dx, j = l,n. (25) 

x 7-i 

It is also convenient to append definition ([M]) with 

fl n = nn - 0-/2X, /xo = ^Sun/2A. (26) 

If the boundary between the Earth's crust and the detector is abrupt, x+ — x~ -C ^ osc , then fi n vanishes. Quite 
analogously, fix vanishes for the abrupt vacuum-Earth boundary. 

By projecting the neutrino state onto the flavor eigenstates, we arrive at the observation probabilities of the 
electron/muon neutrino 

P ^ = { PPe^t) } = 5± iSP^^oVa^^aroM = (27) 
where, for nighttime neutrinos, 

T = r night = isp{a 3 4 ct (x*,x+)a 3 i? d et(x*,x+)e- 2 ^ ^^^a^^ab^^^^a^ 

x g— icr 2 A0n-2 e icrip-n-2 e i<?3Ailin-2 — io^ASi gicipi gi<T 3 Ai/>i i<ri/io g — 2i<7 2 0sun 

x gicri/iog-io-sAi/iigio-ipig-icrgAei e _i<J 3 A^„_i gi<ri/in-i e - i<J 2 Ae„_i e i<Tip„ e -i<T 3 Arp ne ia 2 9~ 1 ^8) 
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while for daytime neutrinos we have 

T = T day = i Sp{a 3 i4 t (xV+)a 3 iW( : EV+)e iCT2 ^ (29) 

In the latter expression, we have omitted the exponential involving /2i , considering the abrupt density change on the 
entry into the Earth. The mixing angle immediately before the detector obviously takes the vacuum value 6q in this 
case. Moreover, the parameter /2o = #Sun/2A is proportional to the inverse radius of the solar core and is thus much 
smaller than the other fij parameters. 

Due to the homogeneity of the detector substance and its smallness compared with the oscillation length, we easily 
find 

R det (x*,x+) = Z+ ot e^ A ^Z- ct = e^-t (1 + {(7 ^ dct)e -^ t + ^2^ (30) 

<T 3 i?dot ( X * ) X n ) °3 -Relet (X*,xt) = 1 - 2i(Tl A^dct Sm 26>det + (<5 2 ot ) , (31) 

where the small parameter 5det is the ratio of the detector width Ldet to the oscillation length £ osc . The quadratic 
remainder terms can obviously be neglected. 



III. DAY-NIGHT ASYMMETRY 



A. Finding the probabilities 



In order to evaluate the probabilities obtained above, let us first make the averaging over the phase Aipi, which 
corresponds to the neutrino path between the creation point xq and the Earth. The region of the neutrino creation 
is extremely large compared with the oscillation length, thus, after this averaging, (cos2a^i) = (sin2A , 0i) = with 
a high accuracy. Using this fact together with Eq. (f3l"|) . after averaging (|29|) we find, up to the terms of the order 
O(fil) = OtfP), 



/gio-3 Ai/>i gicri/io g— 2io-2 0Sungio'iMa g-if3 A^i \ 



( T day) 



COS 2#Sun COS 2flQ W COS 2^Sun, 

cos 29q cos 2#sun- 



(32) 
(33) 



The latter expression coincides with the famous result of Mikheev and Smirnov Q . On the other hand, after averaging 
over Aipi and neglecting 0(r) 2 5 2 ) terms, for the nighttime neutrinos we obtain 



Tnight -> - COS26» S unSp{( 



,2iCT 2 6 



2iffi A^det sin 26» dct )e iCTl p " e i<73 e^ 2 Afl — 1 e iCTl ^ 



ia 3 Aip n - 



x e it73AV'2g-2i(T2Aei e 2i(Ti/ji e — icr3Ai/'2 g— i(T3Ai/^i-i g — i<T2A9„-i gi<Tip„-i g— icrsAVingicriAin | (34) 

where we have made use of the fact that matrices e lai ^ il and e ~ I<J2A6 ' 1 commute up to a negligible term of the order 
O(a0i£i). 

For the calculations which follow, we will use the smallness of the jumps a8i, . . . , A0„_i = 0(f?) and the parameters 
p,i, . . . ,p, n = 0(rj5). Within the linear approximation in the Earth's density parameter 77, this leads to 



Tnight(A#i, . . . , A0„_i;/ii, . . . , /i n ) = cos2#„ COS 



STnjght 



<9T, 



night 



as,u=o 



Partial derivatives with respect to the small parameters are 
<9T nis ht 

Ae,/z=o 



9(A%) 

9T n i g ht 



A6,p,=0 



-icos20 Sun Sp{((7 2 e 2iCT3( >» + 2a 3 AV>dct sin20 dct )e- 2iCT3A ^- } 
2 cos 26*sun {sin 20~ cos 2ai/j„ j — 2Ai/>d G t sin 26 , det sin 2Aifj n j } , 
cos20 S un Sp{(ie 2iCT2e - ai + 2AVdct sin20 det ) e - 2iCT3A ^'} 
2cos26'sun {sin 20 ~ sin2A^ nj + 2A-0dct sin20d e t cos2A?/> n j} . 



(35) 



(36) 



(37) 
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In the above expressions, 

A^„,j = tp{x-) - ip(xj) = \ J u{x)dx = XL nd {l + 0(77)), (38) 

where L n j = x~ — x~ is the distance between the boundary of the jth shell and the detector, measured along the 
neutrino ray. Finally, by substituting the derivatives (j3"6")l and (|3"T1) into Eq. ([33]) and using the fact that sin2A?/Vn = 0, 
we arrive at the final result 

n-1 

Tnight = cos 26> Sun jcos 29~ + 2 sin 29~ ^ (a8j cos 2Aip n j + fij sin 2Atp n j) 

n-1 

-4Ai/>dct sin 26> dot aBj sin 2Aip n<j } , (39) 
3=1 

which is valid up to the terms of the order 0(r]S 2 ) and quadratic terms 0(t] 2 ). The term involving the product of the 
detector width Af/^ct and the oscillating factor ftj cos2Atp n j is of the order 0(r/SS c ict) and is thus omitted. 

The above expression provides a generalization of the main result of paper [25| for the case of nonzero-thickness 
cliffs and detector. It should be stressed, however, that Eq. (f39|) gives poor information on the effects to be measured. 
Indeed, the neutrino experiments last for years, and thus, Eq. (|39j) may only acquire a predictive power after some 
kind of averaging. The averaging procedure should take into account the axial rotation of the Earth (involving the 
integration over the nights), as well as its orbital motion around the Sun. 

B. Averaging the probabilities 

The averaging procedure can be performed analytically, if the oscillation phase incursions Aip n j vary by much more 
than 2ir during the night. Namely, in this case, one can employ the stationary phase technique (see, e.g., [29j]). In the 
case of the beryllium neutrinos (E = 0.862 MeV) traveling through the Earth, the oscillation length is about 30km, 
while the widths of the valleys L n j vary by many hundreds of kilometers, so the oscillation phase variations are indeed 
large enough to use the stationary phase approximation. The only layer, to which this approximation may not apply, 
is the Earth's crust immediately under the detector. This layer is discussed in detail in section IIII CI Appendix I A 41 
and does not interfere with the picture described in this paragraph. 

Let us consider a neutrino traveling through the Earth, which, according to the PREM model [2^, consists of 
a number of concentric spherical shells. The boundary between the valleys Xj corresponds to the point where the 
neutrino crosses one of the interfaces between the Earth's shells; let Tj be the radius of this interface (see Fig. [IJ. 
Further, the distances ■ between the detector and the points where the neutrino enters/leaves the interface with 
radius rj are functions of the 'nadir angle' 6n & [0, n] defined as the angle between the direction to the Sun and the 
nadir in the point of the detector. In terms of the solar elevation angle 9 S [30], the nadir angle is 0n = B s + tt/2. 
The nadir angle, in turn, is a function of the Earth's axial rotation angle r £ [0, 2ir) ('time of day') and the orbital 
motion angle <^ € [0, 2tt) ('season'). The dependence of the distances on the nadir angle is easily found to be 

L n,j = L^jiQz) = r n cos9 N ± ^-r^sin 2 6 N , (40) 
©n ^ arcsin Tj/r n , (41) 

where the upper/lower signs in (1401) correspond to the neutrino entering/leaving the interface rj (see Fig. [IJ. The 
inequality (|4"Tj) ensures that the intersection of the neutrino ray with this interface exists. 

In order to find the night average of the electron/muon neutrino observation probabilities, let us note some properties 
of expressions (|39p and (14TJ1) . First, the number of interfaces crossed by the neutrino is defined via the inequality (|41l) . 
so the number of the valleys and, thus, the number of terms entering the sums in (f39|) are changing during the night. 
The night average of Eq. Q39[) is therefore the sum of the averages of the terms entering this equation and involving 
each of the interfaces j, each average being defined as follows: 

(F(e N )e 2iA ^-) nigllt = / F(e N (T))e 2iA ^ 0N M) . (42) 

J ATnight 
On (r)< arcsin rj / r n 
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FIG. 1: Radial distribution of the electron density N e (r) inside the Earth and the neutrino path through it. The figure 
demonstrates the cross section of the Earth which contains the nadir DO, the center of the Earth O, and the neutrino ray 



Here, F(0n) is some slowly changing function of the nadir angle and Ar n ; g ht is the total duration of the night in 
terms of the Earth's axial rotation angle r, namely, the length of the segment where 0n(t) < n/2 (the Sun is below 
the horizon). 

The duration of the night is, in turn, a function of the season On the equinox, e.g., Ar n ; g ht = 7r, while on the 
winter solstice, Ar n ; g ht — > max. However, the summer nights are just as long as the opposite winter days, so that 



AT„i g ht(^ + 7t) = 2-7T - AT„i g ht(?), 



(43) 
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and the total duration of the nights over all the year is exactly half the year. Therefore, the averaging over the year 
of iV s days should include the division by the total duration of the nights, i.e., irN^. For N q 3> 1, the summation over 
the nights can be replaced by the integration, 



/night, year 



(44) 



and the averaging formula for the terms containing the phase incursion Aip n j takes the form 

dTF(e N (T,0)e 2iA ^ (eN(r < ?)) . 



1 



(F(eN)e 2iA ^'-) night , yC ar = ^2 / <k 



(45) 



6n (t,c) <arcsin Tj jr n 



Now we are able to apply the stationary phase technique to such an integral containing the rapidly oscillating 
exponential. Indeed, let us use the expression, which is valid for smooth functions f(x) and S(x) defined on a segment 
[a, b] containing a single non-degenerate stationary point xq £ (a, b) such that S'(xo) = 0, S"(xq) 7^ [29] 



f(x)e ixsix) dx = 



2tt 



A|S"(z )| 



/(x )exp|iAS'(a; ) + i^ sgnS"'(x )| 



f(y)e l 



i\S'(y) 



+ 0(\- 3 / 2 ), A -> +00. (46) 



The two leading terms come from the stationary point and the boundary, respectively. However, in the application 
to the integral (|4"2"]l . the boundary term vanishes. Indeed, the boundary of the integration domain corresponds to the 
neutrino ray being tangent to the interface with radius rj, hence, d T Arp n j oc d T L n j(Q-^(T)) — > 00, and the boundary 
term is absent. On the other hand, the stationary point is obviously achieved at midnight, when the nadir angle 
0n — > min (the Sun is in its lowest position) , so the integration over the night yields 



7T 

Gn (t) <arcsin rj/r n 



1 



7r\\d*L r , 



r ^(e N )e 2iA ^ Tiw/4 



+ 0((XL n ^-^), (47) 



midnight 



where the two possible signs before i7r/4 correspond to L n j = ■ (see Eq. PO)) 1 ). The principal point here is that 
the second derivative d 2 L n _j at midnight is suppressed for the inner Earth's shells, 



[9 2 L„ J (e N (r))] midnight 
"d£ n j(0 N )" 



a 2 (cose 



d(cos Bn) 



± 



midnight 



dL n ,j(0N) 
d(cos0N) 

sin 2 0m 



midnight 



midnight 



(48) 
(49) 



Now let us use the expression of the nadir ang le ©n via the Earth's axial tilt e = 23.5°, the latitude of the detector 
X € [— 7r/2, 7r/2], and the season q G [0,27r] 30], 

cos 6n (t, <;) = cos x sin ? sin r + cos e cos x cos q cos r + sin e sin x cos <;, (50) 

where ? = corresponds to the winter solstice in the northern hemisphere. The midnight corresponds to the minimum 
value of 8n, which is achieved at r = r m id n ight, 

tan? 



tanr, 



midnight 



(0 



cose 



, cos r midnight (?) = sgn(cos <;) 



cose 



V cos 2 e + tan 2 ? 



Using these expressions, we find the derivative of cos 0n at midnight 

[a 2 (co S e N )] midnight 



cos x cos 2 e + sin 2 ? sin 2 e 



|cos?| v/ C os 2 e + tan 2 q 
Finally, the integral over the night (|47l) takes the form 



f(e N (r)) e 2iA ^- (0N(r)) 



dr 



Gn (t) <arcsin / 



(r 2 /r 2 - sin 2 6 N ) 1/4 
yj 7rAL nj (0N) 



^(e N )e 



2iAj/i„, 3 -=l=iir/4 



(51) 



(52) 



(53) 



midnight 
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On the other hand, the midnight stationary (minimum) values of the nadir angle ON(Tmidnight) vary throughout the 
year (see Eq. (f5"Tj) ). being the smallest on the winter solstice (the darkest midnight) and the largest on the opposite 
summer solstice (the lightest midnight). Therefore, the right side of Eq. (|53|) is still containing a rapidly oscillating 
function of the season and we can perform another isolation of the stationary points, namely, of the two solstices 
S = 0,7r: 

(f (e N )e 2i ^-)„ ig ht, y ea r = / ne N (T,0) e 2iA ^ (eN(T - ?)) ^ 



V 



2n cos e N (T mi dnight(c))| 



©N('?~,<0^ arcs i nT j /r n 
\J r j/ r n- sin2 ©N 

AL hj (0n) 



F(6 N )e 



2iAi/>„j±i(s-l)7r/4 



(54) 

midnight 

2 J sin(e - x)tane, (winter solstice), 

e cose N (T mi dnight(0) = i . . w , , ,. v (55) 

sm(e + x) tane, £ = 7r (summer solstice), 

s = sgn{<9 2 cose N (r mid night(^))} = I ? _ °' (56) 

I "|J-5 S 7t • 

The signs specified in the latter expression are valid in the northern non-tropical latitudes x > e. The detector used 
by the Borexino collaboration is situated in the Gran Sasso laboratory, with 9 = +42.5°, and the prefactor in (|54|) 
which does not depend on j amounts to be 

1 1 1 0.51, winter solstice, 

(57) 



2ny/Af((;)\d% cos 6 N (T m i d night(<0)l 27r ^sin e cos \ | sin(e =j= x) | [0.31, summer solstice 



Therefore, we are left with the following conclusion. The terms entering Eq. (j39]), which contain the oscillating 

functions of the phase incursions 2Aip n _j, are suppressed as O ^ 3 L — = after averaging over the year; 

within the leading approximation, the resulting averages (|54p come from the stationary phase points achieved on the 
winter and the summer solstices. The suppression becomes stronger for the inner Earth's shells. 

It is spectacular that all the terms of the order 0(r]S) in the expression (|3"9")l , including the one corresponding to 
the detector, become 0(r)5 2 rj /r n ) after the averaging. The terms of the order 0(rj), which are proportional to A9j , 
become 0(rjSrj /r n ), respectively. Finally, we are left with the average value 

(T night ) = cos 26 Sun cos 20- + O (nqS^j . (58) 

By substituting this result together with the daytime average value (|3"3"j) into expression (|2T|) for the neutrino obser- 
vation probabilities, we arrive at the day-night asymmetry factor 

A ^ s 2((P e , lght ) - (W) = _^msv,_ . s^ 20 2EV^) + Q t rA 

(-Pe, night) + ("e,day) 1 + ?MSW 008 26*0 Am z \ r n J 

where Tmsw = cos 28q cos 26 , sun = avacasun/wsun defines the observation probabilities for the solar neutrinos due to 
the Mikheev-Smirnov-Wolfenstein effect [2[ and V(x~) is the Wolfenstein potential in the Earth's crust under the 
detector. 



C. The effect of the crust 



The estimation (I5T)1) shown above is substantially based on the piecewise continuous structure of the density profile 
and shows that the asymmetry should depend only on the density of rock immediately under the detector, i.e. in the 
Earth's crust. At the same time, the actual width of the crust is comparable with the oscillation length, and neither 
the valley nor the crust approximation is valid for this layer. Moreover, the phase incursions are small within the 
crust, and so do their time variations, hence, we are unable to use the stationary phase approximation. However, we 
are still able to account for the effect of the crust on the observed day- night asymmetry factor ([5^)1 . relying upon its 
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relatively small density. This feature, together with its bounded thickness, makes it possible to find the closed form 
of the approximate flavor evolution operator for the crust [a;^_ 1; a:^] (see Appendix I A 4p . 

Ro(x-,x+_ 1 ) = cxp{(-ia 2 /3 + ia 1 «)e 2iff ^^- 1 )}+0(7 ? 2 ), (60) 
P + ia = [ Oiyy^^-^-^dy, a, f3 = 0(r]) <E R. (61) 



This result formally repeats the cliff approximation (119[) up to the substitution A0j — ¥ ft, fj,j — ¥ a. Then, exactly the 
same modification arises in (1391) . 



Tnight = cos 26>sun{cos 20~ + 2 sin 26~ ^ (aOj cos 2A^ n ,i + fij sin 2A^ T1 ,j) 

n-l 

-4A^dct sin 29 dct ^ aOj sin 2Ai(> n j j + AT nigh t, (62) 

3=1 

Ar n ig h t = cos26>sun{2sin26>^(/3cos2A'0„ ! „_i + asin2A-0„,„-i) - 4A^dct sm26> do t(/3 sin2A^ n , n _ 1 - a cos 2aVw-i)} 
= 2 cos 20gun sin 29~ J 6(y) cos 2Aip(y) dy ~ 4 cos 26» Su nA^dct sin 26> dct / 6(y) sin 2Atp{y) dy, (63) 

where Atp(y) = ip(x~) — ip(y). The leading 0(77) correction of the crust to the day-night asymmetry factor (p)9"|) reads 

x~ 

A Tmsw sin 2# 2E f ■ 

AA dn = —— / V (y cos2a^ (y )dy. 64 

1 + Imsw cos 20 a Am 2 J 

Again, this result should be subjected to the averaging procedure to acquire the predictive power. However, due to 
the boundedness of the cosine, the correction is easily estimated (both before and after averaging) , 



(A^dn) 



A dn 



Within the Earth's crust, the density change is within 20-30%, hence the change of the Wolfenstein potential 
|(AV) crus t| < 0.25^(2;^). Thus, the non-trivial structure of the crust could produce the deviation of the day-night 
asymmetry factor from expression (I59|) only within 20-30%. The only exception from this rule may occur if the 
detector is placed above the deep region of the ocean (here, both the crust-to-ocean and ocean-to-air jumps of the 
potential are quite large). In this case, the day-night asymmetry factor can be readily found using Eq. f|64[) and may 
become substantially smaller than the one measured by the continental detector. 

IV. DISCUSSION 

From expressions (|59[) and (1641) . one can see that the asymmetry has the order 0{rf) and is determined by the rock 
density in the layer under the detector, i.e. in the Earth's crust. Using the recent data from the SNO and KamLAND 
collaborations 0, namely, tan 2 80 ~ 0.46 and Am 2 w 7.6 x 10~ 5 eV 2 , and the typical electron densities in the 
Earth's crust iV e ( crus t) = 1.3mol/cm 3 and in the solar core N e (s U n) ~ 100 mol/cm 3 [3l|, we arrive at the numeric 
estimation for the day-night asymmetry for solar beryllium-7 neutrinos (E = 0.862 Me V) 

Adn = (-3.8 ± 0.9) x 10" 4 . (66) 

The uncertainty corresponds to the effect of the Earth's crust, which, in principle, can be explicitly evaluated by 
numerical averaging of Eq. (|64l) . Formally, our estimation is in agreement with the results of the Borexino experiment 
on the day-night asymmetry [6], 

A (Borexmo) = ^ ± 12(gtat ) ± 7 ( syst .)) X lO" 3 . (67) 
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The leading corrections not shown in (|66|) come from the averaging over the oscillation phases inside the Earth and 

are suppressed as O (^^r-jf^^j- This suppression is quite strong for dense though deep inner Earth's shells, including 

its core. For the typical number of the valleys n ~ 10, the total contribution of these corrections is several times 
smaller than the effect of the crust specified in (l66l) . 

The energy dependence of the predicted day-night asymmetry factor (JSHJ) within the domain E < 3 MeV is shown 
in Fig. [21 We observe that the asymmetry vanishes for the low-energy neutrinos, which is a result of the applicability 
of the adiabatic approximation for such neutrinos. Indeed, within this approximation, the neutrino flavor observation 
probabilities depend only on the creation and absorbtion points. 

For the neutrinos with energies E > 3 MeV, Eq. (|59|) becomes only qualitative. For instance, the energies of the 
boron neutrinos, which were considered, e.g., in the paper [25[, reach 10 MeV, while the oscillation lengths are as large 
as 300 km. For such energies, both the relative errors of the valley solutions (fl~8")l (multiplied by their number n ~ 10) 
and the averaging errors may, in principle, become quite large. 

One should also mention here that, although the day- night asymmetry factor is a quantity which is directly measured 
in the neutrino experiments, both its sign and magnitude depend not only on the neutrino regeneration effect inside 
the Earth, but also on the value of the Wolfenstein potential in the solar core, where the neutrino is created. This fact 
is easily seen in (|59[) . where the solar effect manifests itself via the quantity Imswj which depends on the neutrino 
energy E. At the energy E ~ 2.3 MeV, which corresponds to the Mikheev-Smirnov resonance, we have as U n = and 
Tmsw = 0, and, as a consequence, the day- night asymmetry factor (|59[) changes sign (see Fig. [2]). Thus, it is the 
resonance inside the Sun which makes the asymmetry, observed on the Earth, vanish. 

AdnXlO 3 




FIG. 2: The predicted day-night asymmetry factor Adn as a function of the neutrino energy E plotted for Am 2 = 7.6 x 10" 5 eV 2 
and tan 2 8 = 0.46 

Finally, one may observe that beryllium neutrinos (E = 0.862 MeV) are indeed quite useful for the study of the 
matter effects in the neutrino oscillations, since they correspond to almost maximum negative asymmetry factor. 
Moreover, solar beryllium neutrinos are highly monochromatic (in contrast, e.g., to the boron neutrinos) and their 
flux is considerably larger [3l[ . Thus the Borexino experiment is a prospective attempt to detect the day- night flavor 
asymmetry. However, we are able to conclude that the latter effect needs a 10-20 times better experimental resolution 
to be distinguished at a considerable confidence level. This conclusion does not depend on the densities of the inner 
Earth's layers and on the substance and the dimensions of the detector, as far as they arc considerably smaller than 
the oscillation length. 

As a demonstration of our conclusions about the time average of the day-night asymmetry factor, we present the 
results of its numerical averaging, together with the theoretical curve (|5fJ)) . in Fig. [3] The two dotted curves correspond 
to the neutrino detectors placed in Gran Sasso (x = +42.5°) and on the Northern tropic (\ = +£ = +23.5°). The latter 
curve demonstrates specific oscillations which come from the amplified contribution of the winter solstice stationary 
point to the year average of the asymmetry factor (see Eq. ([54)0 . Moreover, the approximate smoothness of the both 
curves indicates the (approximate) stability of the day-night asymmetry factor with respect to slight variations of the 
parameters of the PREM model, namely, the radii of the Earth's shells and the density jumps. The stability becomes 
stronger for low-energy neutrinos, as well as for the detectors operating far from the tropic. 

It is also interesting to study the effect of the local Earth's crust inhomogeneities under the detector on the observed 
day-night asymmetry. Such inhomogeneities could be associated, for instance, with the oil-bearing horizons. Let the 
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FIG. 3: Numerically averaged day-night asymmetry factor versus the predicted one for different neutrino energies E. The 
curves with circle and triangle fit points correspond to the neutrino detectors placed in Gran Sasso (x = +42.5°) and on the 
Northern tropic (x = +£ = +23.5°), respectively. The solid curve corresponds to the theoretical estimation (|59[) of the average 
asymmetry 



inhomogeneity be described by the variation SNe (x) of the electron density over the smooth profile N e (x), 

N e (x) = N e (x) + 5Njp(x), SNP(x) = for x g [x h Xi + Sxi], (68) 

where the inhomogeneity size 8x\ <C ^osc- Then the contribution of this inhomogeneity to the asymmetry factor is 
given by Eq. 



SA 



(0 

dn 



1 



N, 



(crust) 



5 'N® '(y) cos 2 Aip(y)dy = 



1 



Xi-\-5x{ 



N. 



e(crust) 



SA 



A, 



< 



\SNP\ 2nSxi 

-^e(crust) ^osc 



2?rL( i ) 



sm ■ 



5N?\y)s^2^{y) 2 -^^, (69) 



(70) 



where the depth of the inhomogeneity under the detector. Therefore, exploration of the Earth's 

crust based on the neutrino oscillations would require an extreme improvement of current measurement techniques 
[32| . and its prospects seem obscure in the nearest future. 
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Appendix A: Approximate solutions for the evolution operator 



In this section, we derive the approximate solutions for the flavor evolution operator Rq(x, X q) (see Eqs. © and 
(fT4")) ) for the valleys and cliffs. First, due to the relatively small Earth's density, which manifests itself as a small 
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parameter 77 < 10 2 , we resort to the linear approximation in 77, namely, take the two leading terms of the Dyson 
series 



R (x,x ) = Texp<( -ia 2 j 9(y)e 2[a3 ^ v) dy \ = 1 - icr 2 J 0(y)e %a3 ^dy + 0{r] 2 ). 

Xq ) X 



(Al) 



1. Valleys 

In the valleys, we have a smooth and bounded function 9 (y) and a rapidly oscillating matrix exponential e 2 " 73 ^^^ = 
cos2ip(y) + ia 3 sin2^(y). Then, applying the double integration by parts, we arrive at 



X 

I 



9(y)e 2ic " 3,p( - v) dy = ~ V y) e (v) c 2ia 3 ^( y ) 
2ia 3 \uj(y) 



V 2 y 9{y) ■ e 2i(73 ^dy 







v vf(y) = it. 



1 



dy \2ia 3 Xuj(y) 



(A2) 



(A3) 



Since the total variation of the effective mixing angle 9(x) is bounded within the valley and has the order 0(77), all 
the gradients in the above expressions are suppressed by powers of the small parameter S = £ osc /XL = Tt/\L, where 
L is the width of the valley and £ osc is the oscillation length. Namely, 



0(a9/L) = 0( V /L), 
0(a9/L 2 ) = 0( V /L 2 ), 



uj = y/a 2 + b 2 = l+ 0(rj), 

= Hp = OWL). 



Using these estimations, one can readily show that 

l('9 ibf 



0(r,S/L), 



V 2 y 9(y) = Oi^/L), 
then the matrix norm of the remainder integral in Eq. (|A2 



V 2 9\y) ■ e 2ia ^ iv) dy 



< / \\Vp(y)\\dy = 0( V 5'). 



(A4) 
(A5) 
(A6) 
(A7) 



(A8) 
(A9) 

(A10) 



On the other hand, the first term in the right side of Eq. (|A2[) consists of two parts, the one proportional to T> y 9(y), 
which is of the order 0(r]8/\L) — 0(r)5 2 ), and the other proportional to 9(y), which is 0(r)5). Then, summarizing 
the estimations made, we conclude that 



; 2i<r 3 ^(y) i _ 



2ia 3 \uj{y) 



= 2ier 3 V>(y) 



0( V 5 2 ). 



(All) 



Moreover, here, due to (|A6j) . we can safely replace uj(y) by unity, and then the approximate solution of the evolution 
equation in the valley takes the form 



R (x,x ) = 1- ^ (9{x)e 2i ° 3 ^ x) - 9(x Q )e 2i(73 ^ Xo ^ + 0(rjS 2 ). 
Finally, by neglecting the terms of the order 0{rj 2 8 2 ) 1 we may write 

R (x, xo) = exp (9{x)e 2i ° 3 ^ - 9(x )e 2i ° 3 ^) J + 0( V S 2 ). 

This expression can also be derived using mathematically strict stationary phase technique [2£ 



(A12) 



(A13) 
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2. Cliffs 

On the cliffs, the change of the effective mixing angle is still of the order 0(r]), while the spatial scale of this change 
is quite small. On the other hand, the change of the phase of oscillations Aip <C 2tt, so that we can expand the 
exponential in the right side of Eq. (|A1[) in the local phase incursion 

v 

V>(y) - $(x ) = xj u(z)dz = X(y - x ) + 0(r)6), (A14) 

x 

e 2ia 3 ^ y ) _ e 2i^(»o)( 1 + 2i <73 A(j / -a;o)) + 0(8 2 ) + 0(r]S). (A15) 
Now, using the smallness of the total variation of the mixing angle on the cliff, we obtain 

X 

Q(y) e 2w1>(v)dy = e 2i<T 3 V(xo) f ^ + 2ia 3 X(y - x ))6(y)dy + 0{r,5 2 ). (A16) 



This finally leads to the cliff approximation of the evolution operator 



Ro{x,x ) = 1 + {-ia 2 A0 + ia 1 ii)e 2ia ^ {x ° ) + 0{ V S 2 ) 

= cxp{(-ia 2 Ae + ia 1 ^)e 2i ' 7 ^ < - x ^} + 0(t 1 6 2 ) + 0(t] 2 ), (A17) 



where a8 = 6(x) — 0(xq) = 0{rj) and 



H = 2X(y- x o )0(y)dy = 0( V S). (A18) 



XO 



3. Valley + cliff 

It is also useful to consider a valley [Xj,xJ +1 ] following the cliff \xj,x^]. Using expressions (|A13|) and (|A17j) for 
the evolution operators on these segments, within the linear accuracy in 77, we obtain 

R (xj +1 ,xJ) = Ro(xJ +1 ,xj)Ro(xj,xj) 

= exp {e{xJ +l )e 2i "^+0 - 6(x+)e 2i ^ (x t^ + (-kr 2 A0j + ia^e 2 ™*^ J + 0{ V S 2 ) 

= exp j-|p(z- + J e 2iCT3 ^ + i) + (~ia 2 A6 3 ; + kn^) e 2 ^^ J + 0( V S 2 ), (A19) 

where A0j = O(x^) — 6(xj), fij = [ij + 8(x^)/2X, and \ij is defined analogously to (|A18|) . Now, by substituting the 
above result into representation (j9)) for the evolution operation, we arrive at 



R(xJ +1 ,xj) = e iCT2e ^ +1 )e iCT ^ ( ^ +1 )i? (x7 +1 ,x+)i? (4,a;-)e- iCT ^ ( ^)e- iCT2e ^) 



J<? 2 0(x j + 1 ) 



-i^e(xj) + ^ 2 j 

= e ia 2 e(x7 + J e -ia 1 e(x7 + J/2A e ia 3 A^ +le ia 1 /i Je -ia 3 Ae 3e -i<r 2 0(a : 7) + (r]8 2 ) + 0(jf), (A20) 

where A^j+i = VK^J+i) — ^{ X J)- This expression, in turn, can be generalized to a sequence of n valleys (see Eq. (|2"T)| ). 



4. The Earth's crust 



As one could see from the main flow of the paper, for the neutrinos with energies around 1 MeV, the Earth's 
crust cannot be considered either a valley or a cliff, since its width is comparable with the oscillation length, and the 



15 



parameter S is of the order of unity. Here, however, another approximation is useful, which takes into account the 
small density of the crust, as well as its bounded thickness. In terms of the parameters rj, and 5, we have 

?7crust = 2EV "^ st ~2x 1CT 3 , 5 crust = -7— < 5 (beryllium neutrinos). (A21) 

A 771 t osc 

Then, using the expansion (|A1[) . we find the immediate expression for the evolution operator for the crust: 

R Q (x, x ) = (1 + ioia - ia 2 f3)e 2i,73 ^ Xo} + 0(r] 2 ) = exp{(-ia 2 /3 + ia 1 a)e 2i ' J ^ {xa) } + 0(t? 2 ), (A22) 
where real numbers a, (3 = 0(rj) are defined by the expression 

X 

P + ia= [ e(y)e 2i ^-^ X0 ^dy. (A23) 



Note that the form of this approximation (|A22|) coincides with the cliff approximation (|A17|) . up to the coefficient 
substitution aO — > (3, \i — > a. In particular, the cliff approximation is restored in the 5 — > limit. 
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